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Abstract 



Based on the correspondence of projection to an eigenvalue problem via the 
underlying group structure, we generalize the projection algorithm proposed 
by Morrison et al. This algorithm is eligible for large-scale computations, 
because we can avoid the accumulation of errors. As an example, we apply it 
to the spurious cm. motion in a multi-ftu; shell-model calculation. In com- 
parison with Lawson's method, the present method is exact and nevertheless 
does not need longer CPU time. It is pointed out that the exact treatment of 
the spurious cm. motion is crucial to some physical quantities. 

PACS numbers: 02.70.-c, 21.60.Cs 
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In the numerical treatment of a Hamiltonian in a large space, we often need an efficient 
algorithm for projection, in connection to conservation laws or spurious motions. For in- 
stance, the J (angular-momentum) projection accelerates shell-model calculations in many 
cases. A more serious problem is the spurious center-of-mass (cm.) motion: since we some- 
times use a Hamiltonian which may excite the spurious motion, it is important to project 
it out of basis vectors. In this report we extensively develop a projection algorithm, which 
was first proposed in Ref. |I| for the J projection, and apply it to the cm. motion problem 
within the regime of multi-^o; shell-model calculations. 

Suppose that, in an iV-dimensional (N is finite) Hilbert space V, we have relevant and 
irrelevant components. The projection operator P picks up only the relevant components, 
eliminating the irrelevant ones. We now introduce a certain hermitian operator Hp, which 
satisfies HpV C V and has degenerate eigenvalues within V. We denote the eigenvalues by 
Ho, Hi, ■ • ■ , Hn-i (n < N). Eigenvectors associated with each eigenvalue Hk span a subspace 
Wk (k — 0, 1, • • - n — 1). The full space V is equal to the direct sum of Ws; 

n-1 

V = ®W k . (1) 

fe=0 

Any vector \v) E V is expanded as 

n-1 

H = E^ (fc) >, (2) 

k=0 

where \v^ k >) E Wk (i.e. Hp\v^) = fik\ v ))■ It is postulated that the relevant components 
belong exclusively to Wo, which is attached to the eigenvalue /io; PV = Wo and P\v) oc \v^). 
(We can straightforwardly extend our discussion to the case that the relevant subspace is 
a sum of a few W's.) The projection can then be redefined as an eigenvalue problem with 
respect to Hp, as will be illustrated below. 

It is crucial to find such an operator Tip with a simple form. The projection is normally 
linked to a certain group structure. As is well-known, conservation corresponds to a sym- 
metry, and therefore to a group. A spurious motion occurs when we ignore a symmetry to 
be satisfied. This symmetry leaves its trace in a group structure in the enlarged space. In 
both cases, the relevant components are characterized by a specific representation of the 
corresponding group. Then a Casimir operator of the group can be adopted as Hp; for 
instance, Hp = J 2 in the J projection. The condition HpV C V is required in order that 
the group structure holds within V. 

For a given vector \v) E V, \v^) is obtained in the following manner, by applying the 
Lanczos diagonalization algorithm with respect to Hp. Having a set of basis vectors T n i = 
{\v),Hp\v), ■ • • , (Hp) n with the Gram-Schmidt orthogonalization, Hp is represented 

by a tridiagonal matrix. We denote the component of this matrix by (Hp)ij. At a 
certain n'(< n), (Hp) n \v) is no longer linearly independent of Y n i because of the degenerate 
eigenvalues. Then \v ^) is obtained by solving a system of coupled linear equations 

£(Wp)iizf = Hkxf ) (i = l,...n'), (3) 
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where x\ yields the expansion coefficient of \v^} for the Lanczos bases and = Ck of 
Eq.(||). Since Hp is tridiagonal now, the sum over j on the left-handed side of Eq.(0) runs 
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from i — 1 to i + 1 at most. What we need is the eigenvector \v^). The probability (v\P\v) 
is given by |£i°' > | 2 . It has been known that solving Eq.(||D from i — n' to 1 with a small initial 
value of $ n / is advantageous, in order to keep rounding errors small @. 

Because of the group structure, Hp has highly degenerate eigenvalues; n <C N . For the 
J projection in shell-model calculations, n is the number of possible J values in the space 
V, being several tens even for the M-scheme shell- model space with N ~ 10 6 [H]. In the 
above algorithm, we can keep the number of bases small (n' < n). As far as we know Ho in 
advance, the inverse iteration method is immediately applicable to obtain |i>(°)). However, 
solving Eq.(^|) from i — n' to 1 is more efficient, because we reach \v^) with a small number 
of operations. The Lanczos method has an advantage for these reasons, although n is not 
quite large. 

It is noted that, by using Hp, we can explicitly construct the projection operator P as 

p = g(^zB). (4) 



k=l 



HO — 



Each factor (Hp — fik)/(no — Hk) eliminates the irrelevant component \v ^) G Wk with k ^ 0, 
leaving |i> < -°' ) ) at last. Despite its mathematical exactitude, this procedure is not numerically 
efficient. First of all, we always need an n-fold operation of Hp, in contrast to n'(< n) in the 
present algorithm. Moreover, in the process eliminating \v ( k >), other irrelevant components 
may be enhanced, causing serious numerical errors. We take the J = projection as 
an example, with Hp = J 2 and = + 1). We operate (Hp — fik)/{po — Hk) on 
\v) from k = 1 to n — 1. The irrelevant component |t»W) is removed at the first step. 
However, if a small amount of the W„_i component is contained in \v), this component is 
enhanced by the factor [n(n — 1) — 2]/2 ~ 0(n 2 ). After several steps, the W n -i component 
is repeatedly enhanced and eventually dominates over the relevant component, giving rise 
to a loss of precision for the relevant component in practical floating-point computations. 
It is numerically advantageous to enhance the relevant component, which corresponds to 
a single eigenvalue /j,q, rather than to remove many irrelevant components. In the present 
algorithm, the relevant component is enhanced because of the order of solving Eq.@. 

Morrison et al. have pointed out the equivalence between projection and an eigenvalue 
problem for some cases, particularly for the J projection [0. They proposed the Lanczos 
algorithm with Hp = J 2 . However, it has not been discussed how to obtain Hp in general 
cases, since one has not fully realized the underlying group structure in the correspondence 
of the projection to an eigenvalue problem. Furthermore, only the case of [H, Hp] = was 
discussed in Ref. UjJ. Generalizing their argument, we have now developed the projection 
algorithm systematically. 

Since we apply the Lanczos method, Hp as well as H do not have to be represented in 
a matrix form (|]. The present algorithm may be extended even to the case that n (and 
therefore N) is infinite or the eigenvalues of Hp are continuous. In such cases, however, the 
degree of convergence with finite n' strongly depends on the vector \v). We have restricted 
ourselves to the finite iV cases in this article, by which the convergence for n' is assured. 

We next demonstrate the present projection method for the spurious cm. motion in 
a mvlti-huj shell-model calculation. Shell-model calculations in mvlti-Tiu spaces are under 
current interest, in relation to a microscopic description of nuclei and to a research on 
light unstable nuclei 0. Because the origin of the coordinate space is fixed in the shell 
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model, the cm. motion, which is irrelevant to the nuclear structure itself, is included in 
the model space. With the harmonic-oscillator (h.o.) single-particle wavefunctions, the cm. 
motion is also expanded in terms of the h.o. representation, which forms a U(3) group. It 
is desired to constrain the cm. mode to the Os orbit. For this purpose, we take 

1 - AMu 2 n7 3^ 
~ 2AM 2 _ 2 ' < 5 » 

where R and P denotes cm. coordinate and momentum. This is indeed proportional to a 
Casimir operator (linear Casimir operator) of the U(3) group, and is expressed with up to 
two-body operators within the shell-model space, 
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Mu? 2 





Pj + Muj 2 Yi ■ r. 



(6) 



The Os mode corresponds to the eigenvalue /io = 0. 

In order for the cm. motion to be well-defined inside the model space V, HpV C V 
should be satisfied. Otherwise the model space should be extended. We consider a model 
space with the hu-type truncation; V is a direct sum of subspaces characterized by the 
number of oscillator quanta m, V = (& m V( m \ Each subspace V^ 71 ' is expanded by direct 
products of the relative motion and the cm. motion, both of which have specific numbers of 
oscillator quanta, = ® k [V™l k <g> W£ m -]. A resummation as W k = ©^Vf 1 ® W^ m '] leads 
to the expression of Eq.(|l]). The space W k is characterized by the symmetric representation 
[k] of the cm. U(3) group and has = kTiu). Obviously HpV C V is satisfied. All the 
components outside Wq are regarded as spurious components. The corresponding projection 
operator P of Eq.([|) has been discussed in Ref. 0. Unlike the J projection, we do not 
necessarily demand that the shell- model Hamiltonian H should commute with Hp ||, under 
the presence of a core or inactive particles. 

The present projection is incorporated into the shell-model diagonalization. The Lanczos 
method is usually employed in diagonalizing a shell-model Hamiltonian The diago- 

nalization of H within Wq is equivalent to the diagonalization of PHP in V, except that 
the irrelevant states have null eigenvalues for PHP. Note that the explicit form of PHP is 
complicated. Given an initial vector \vi), we apply the present projection and obtain |wj ^}. 
The next Lanczos basis \v% ) is generated from PH\vi }. Multiplication by H is performed 
just as in the usual shell-model calculation. The projection on \v 2 ) = H\vi) yields l^ ^), 
from which \v% ) is obtained by orthogonalizing to {v^). The Lanczos bases of PHP are 
generated one after another in this manner, as is summarized in Fig. [I]. The computation 
proceeds according to the arrows. Each Lanczos diagonalization is marked on the boxes. At 
left the main process of the diagonalization of PHP is indicated, while the projection by 
diagonalizing Tip is shown at right. 

When the single-particle energies consist only of kinetic energy and the effective interac- 
tion depends only on relative coordinates with the h.o. bases, the shell-model Hamiltonian 



can be made commutable with Hp |ID| . Even in such cases with [H, Hp] = (i.e. [H, P] = 0), 
the present projection is useful to cut off numerical errors. It may be applied in the Lanczos 
iterations of if at a certain interval, as well as to the initial vector. 

For the spurious cm. motion problem, the so-called Lawson's method has conventionally 



been used [i,p,10|. The original Hamiltonian is amended as 
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H' = H + \n P (A > 0) 



(7) 



The Tip term shifts up the spurious components and reduces their admixture in low-lying 
states. The quality of this prescription for large A can be assessed by perturbation theory. 
By regarding XTip as an unperturbed Hamiltonian and if as a perturbation, an eigenvector 
of H' is approximated by 

n— 1 i 

\v x ) - |t,(°>) - Evk\v (k) ) ; Vk = —( v W\H\v^) ~ 0(A -1 ) . (8) 

As far as [H, Tip] ^ 0, r\ k does not vanish. Notice that rjx is the largest in general, because 
of the energy denominator. The corresponding energy is 

n-1 

E x = (v x \H'\v x ) - E® - £ Vh (vV>\H\vM) + 0(A~ 2 ) , (9) 

fc=i 

where f^ -* = (t> ^\H\v ( -°* ) ). Thus the influence of the spurious modes on the wavefunctions 
and energies is proportional to A -1 , not diminishing rapidly as A increases. On the other 
hand, as we set larger A, the convergence with respect to H' becomes slower, because the 
Tip term dominates and fine splitting due to H is not easily attained. Moreover, even if 
the influence of the spurious motion is small for energies, it is not necessarily negligible for 
transition properties, particularly when there is some cancellation. 

We have carried out a numerical test for the 10 Be nucleus. In the lowest configuration we 
have six nucleons in the p-shell, on top of the 4 He core. We take the model space V spanned 
by up to the 3hu excitations (n = 4). For H, the single-particle energies are determined from 
the Woods-Saxon potential with the parameters of Ref. ||11|| . The Millener-Kurath [1^] and 



USD [O] two-body interactions are taken for the (s-p-sd)-she\\, while the M3Y NN effective 



interaction |14| is used for the rest. This Hamiltonian is adopted just for a numerical test, 



putting aside its physical appropriateness. As a matter of fact, the negative-parity levels 
have too high excitation energies. 

As in Fig. [1], the subroutine for the cm. projection is incorporated into the shell model 
code VECSSE 0,[[5|. The present projection technique is compared with Lawson's method. 
For the A parameter in Lawson's method, we take (i.e. no care of cm. motion), 1, 5 and 20, 
with hu = 45A- 1 / 3 -25A- 2 / 3 = 15.5 MeV. The CPU time on HITAC S3800/480 is presented 
in Table [IJ Both the positive- and negative-parity states are computed up to J = 4, and 
5-7 lowest eigenstates are computed for each J p , with the convergence criterion of 0.1 keV. 
While the present method gives an exact projection, it does not make the CPU time so 
long. In Lawson's prescription, it takes a longer CPU time as A increases, because more 
iterations (i.e. more Lanczos bases) are necessary to fulfill the convergence criterion. While 
for A = 1 the iteration number is almost the same as that in the present projection, it costs 
almost twice as many iterations for A = 20, resulting in longer CPU time than the exact 
projection. The calculated energies are presented in Table y. The energy eigenvalues of 
Lawson's method do not rapidly converge for increasing A; the 100 keV discrepancy remains 
even with A = 20. The excitation energies, on the other hand, are rather stable. For A = 0, 
the spurious component dominates the lowest 1~ state with E x = 13.998 MeV, while the 
third lowest 1~ state consists mainly of the Wq component, whose energy is presented in 
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Table [T| As a measure of admixture of the spurious components, the expectation value of 
Tip is also shown for 4f. 

We turn to the transition properties. Some of the calculated results are presented in 



Table [II. The El transition strengths are rather susceptible to the contamination of the 



spurious cm. motion, since the bare El operator T(E1) = eJ 3/47rSie P r i has the isoscalar 
part proportional to R. Here i"i denotes the coordinate of the i-th nucleon and the sum runs 



over all the protons. A simple modification of T'(E1) = eJ3/4:irJ2i£ P ( r i ~ R-) reduces the 



influence of the spurious motion ||16|| . However, we do not find apparent improvement for 
Lawson's method in Table [II. On the other hand, the present method yields equal values 
between T and T', since the cm. excitation is entirely removed. 

A similar treatment for the E2 transition to T'(E1) requires two-body operators, and 
hence is not popular. We here assume the usual one-body E2 operator with the bare charges. 
Compared with the exact values, Lawson's method works well in most cases. However, the 
I2 — > 24 transition, whose strength is not large, shows a significant discrepancy. Even when 
excitation energies seem good for a physical discussion, the influence of the spurious motion 
may be serious for transition strengths, particularly for relatively small ones. 

In summary, we have developed a projection algorithm, extending the method of Ref. 
|Ij]. The projection is redefined as an eigenvalue problem of a Casimir operator of the 
underlying group. On this ground, the Lanczos diagonalization technique can be utilized also 
for projection. This algorithm is efficient particularly in large-scale computations, because 
numerical errors do not accumulate severely. We have applied it to eliminate the spurious 
cm. motion in a multi-huj shell-model calculation. Whereas the present method is exact, it 
does not need longer CPU time than Lawson's method. The present algorithm is useful, or 
even crucial, to evaluate physical quantities sensitive to admixture of the cm. excitation. 

We thank A. Gelberg for careful reading the manuscript. The Millener-Kurath and USD 
interactions registered in the OXBASH program package [17] have been used, with the 
permission of B. A. Brown. 
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TABLES 



TABLE I. CPU time (sec) for diagonalization of H, with the cm. motion handled by the 
present algorithm or Lawson's method. 



Projection 




Lawson's method 






A = 


A = 1 A = 5 


A = 20 



586 523 535 560 629 



TABLE II. Eigenenergy of the ground state and excitation energies (MeV). {Hp) value is also 
shown for 4+ in the parenthesis. 



Exact A = A = 1 A = 5 A = 20 

E(0T) -375.845 -376.691 -376.453 -376.142 -375.950 

E x (2f) 4.246 4.200 4.212 4.229 4.240 

E x (4f) 17.922 17.617 17.707 17.821 17.888 

{{Hp)/flio) (0.0) (1.3 X10" 1 ) (1.5X10- 2 ) (3.3 x ltr 3 ) (3.9 xlO- 4 ) 

E x (l{) 19.266 19.125 19.166 19.263 19.269 

E X (S^) 19.409 18.234 19.476 19.445 19.422 



TABLE III. B{E1) (e 2 fm 2 ) and B(E2) (e 2 fm 4 ). 







Exact 


A = 


A = 1 


A = 5 


A = 20 




1 2 ) by T 


0.901 


0.292 


0.953 


0.888 


0.898 




hyT' 


0.901 


0.498 


0.690 


0.878 


0.896 


B(E2;0+ 


-2+) 


13.167 


12.905 


12.949 


13.045 


13.121 


B(E2; 17/ 


-3D 


2.059 


1.807 


2.047 


2.059 


2.059 


B(E2; 1 2 " 


-27) 


0.011 


0.413 


0.265 


0.002 


0.004 
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FIGURES 
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FIG. 1. Illustration of the present projection algorithm incorporated into the Lanczos diago- 
nalization of H. The basis vectors are generated so as to stay inside Wq. 



9 



